None
import tifffile as tf
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
from scipy import spatial
import pandas as pd
import seaborn as sns
import PIL
from math import sqrt
import matplotlib.patheffects as PathEffects
import random
from tqdm import tqdm
from scipy.optimize import minimize
from PIL import Image
from os import listdir
from os.path import isfile, join
import warnings
from helper_functions import *
#this alters the amount of width taken up but notebook chunks
from IPython.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
well_ID = 12
lomag_tiles_dir = '/path/to/your/raw/low_magnification/tifs'
himag_tiles_dir = '/path/to/your/raw/high_magnification/tifs'
#replace these with your paths at runtime
lomag_tiles_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/input_files/raw/genotyping_tifs/cycle_0/'
himag_tiles_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/input_files/raw/phenotyping_tifs/Mito_wells/'
lomag_df = RetrieveCoordinates_in_df(lomag_tiles_dir, well_ID)
himag_df = RetrieveCoordinates_in_df(himag_tiles_dir, well_ID)
lomag_df.head()
#calculate radius and centroid of higher magnification tiles
_well_radius, _well_centroid = ApproximateRadiusAndCentroid(himag_df)
well_fiducial_tiles = GenerateFiducials_df(lomag_df, himag_df, _well_radius, _well_centroid)
well_fiducial_tiles.head()
#both dataframes are generated from OPS outputs from processing pipelines documented in this GitHub repo: feldman4/OpticalPooledScreens
#replace these with your output at runtime
df_lomag = pd.read_csv('OPS_genotyping_output.csv', index_col=0) #readcalling output from OPS pipeline
df_himag = pd.read_csv('OPS_phenotyping_output.csv') #output from OPS pipeline
df_lomag.sample(3)
df_himag.sample(3)
#outputs from standard OPS processing of raw images
lomag_processed_images_dir = '/path/to/your/OPS_processed/low_magnification/tifs'
himag_processed_images_dir = '/path/to/your/OPS_processed/high_magnification/tifs'
#replace these with your paths at runtime
lomag_processed_images_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/Readcalling/Mito_process_05072023_thresh90_v1/'
himag_processed_images_dir = '/home/sansburys/shalemlab/v2_OpticalSequencing/data/VL1_JFX549_3May2023/Phenotyping/Mito_process_05072023_pheno_v1/'
#instantiate empty lists to fill with local coordinate information
cellpairs_list = []
tilepairs_list = []
#for each pair of tiles in the well_fiducial_tiles dataframe, print out images and identify the number associated with the same cell at each magnification
PrintImagesWithAllNucCentroids(tile_lomag=157, tile_himag=2385, lomag_image_dir=lomag_processed_images_dir, himag_image_dir=himag_processed_images_dir,
cells_lomag_df=df_lomag, cells_himag_df=df_himag,
well=well_ID, suppress_labels=False, sample_lomag=0.8, _vmax=2500)
#after investigation of above images, supply manually identified fiducial cell identities
#cell_IDs
cellpair_lomag_himag = [174,12] #of form [cell_ID_in_lomag_image, cell_ID_in_himag_image]
cellpairs_list.append(cellpair_lomag_himag)
#tile_IDs - these will typically match the tileIDs provided in the dataframe, but on occaision, you may choose to search neighboring tiles if a shared cell can't be identified
tilepair_lomag_himag = [157,2385]
tilepairs_list.append(tilepair_lomag_himag)
#repeat this for every line in the well_fiducial_tiles dataframe
#after repeating this for 8-12 fiducials total:
cellpairs_list = [[1683, 36], [160, 8], [461, 19], [863, 22], [2587, 97], [1911, 92], [17, 16], [129, 20], [2063, 62], [45, 8], [675, 22], [655, 12]]
tilepairs_list = [[8, 64], [157, 2385], [89, 1227], [76, 1171], [32, 298], [23, 334], [133, 2115], [142, 2148], [40, 499], [125, 1947], [86, 1238], [79, 1182]]
fiducial_cells = pd.DataFrame({
'cell_lomag': [pair[0] for pair in cellpairs_list],
'cell_himag': [pair[1] for pair in cellpairs_list],
'tile_lomag': [pair[0] for pair in tilepairs_list],
'tile_himag': [pair[1] for pair in tilepairs_list],
'well': well_ID})
#look up tilewise coordinates for each of these cells to later place on global coordinate scale
points_hi, points_lo = GetTilewiseLocations(df=fiducial_cells, himag_df=df_himag, lomag_df=df_lomag, well=well_ID)
points_hi, points_lo
# Get image height and width
example_tif_lomag = os.path.join(lomag_tiles_dir, random.choice([f for f in os.listdir(lomag_tiles_dir) if f.endswith('.tif')]))
example_tif_himag = os.path.join(himag_tiles_dir, random.choice([f for f in os.listdir(himag_tiles_dir) if f.endswith('.tif')]))
h_lomag, w_lomag = Load_TIF(example_tif_lomag).shape
h_himag, w_himag = Load_TIF(example_tif_himag).shape
# Generate a starting grid to represent how tiles relate to each other in 2d space.
# This will depend on the manner in which the stage moves as images are collected, and will be investigated later and corrected if necessary
untested_grid_lomag, lomag_coord_df = GenerateTileGrid(path_name=lomag_tiles_dir, well=12)
untested_grid_himag, himag_coord_df = GenerateTileGrid(path_name=himag_tiles_dir, well=12)
print(untested_grid_lomag)
#let's test arrangment by picking a set of four tiles that are adjacent in the grid, and determining which transformation (if any) is required to align them
input_dict = {'top_left':81,'top_right':82,'bottom_left':98,'bottom_right':97}
TestTileGrid(input_dict, lomag_tiles_dir, lomag_coord_df, zoom=True, zoom_factor=0.7) #it can be helpful to zoom in when checking lomag images
#confirming same is true for high magnification images
input_dict = {'top_left':439,'top_right':440,'bottom_left':512,'bottom_right':511}
TestTileGrid(input_dict, himag_tiles_dir, himag_coord_df, zoom=False, zoom_factor=0)
#update tile grid to reflect most seamless orientation
arranged_grid_lomag = SetTileConfiguration('flipud', untested_grid_lomag)
arranged_grid_himag = SetTileConfiguration('flipud', untested_grid_himag)
#Use updated tile arrangement grids, pre-calculated fiducial points at each magnification, and image dimensions to calculate mapping parameters
P_lomag = Local_to_Global(points_lo, arranged_grid_lomag, [h_lomag, w_lomag])
P_himag = Local_to_Global(points_hi, arranged_grid_himag, [h_himag, w_himag])
mapping_parameters = Fit_By_Points(P_lomag, P_himag, verbose=True)
df_all_reads = df_lomag
subset_cells = df_all_reads[(df_all_reads["design"]=="ENST00000309318.8_1_3__11") & (df_all_reads["well"]==well_ID)& (df_all_reads["cell_barcode_count_1"]<1)].sample(10)
lomag_subset = np.array(subset_cells[['tile','i_nucleus','j_nucleus']]).astype(int)
#calculating global coordinates from lomag image
calculated_P_lomag_subset = Local_to_Global(lomag_subset, arranged_grid_lomag, [h_lomag, w_lomag])
#using mapping parameters calculated above to calculate global higher mag coordinates
calculated_P_himag_subset = model_TRS(calculated_P_lomag_subset, mapping_parameters, angle='degree')
#use global higher mag coordinates to translate to local higher mag coordinates
reverse_calc_p_himag_subset = Global_to_Local(calculated_P_himag_subset, arranged_grid_himag, [h_himag, w_himag])
#use local lomag coordinates and calculated local higher mag coordinates to look up the same cell at each magnification
#inspect visually that red dot lands on same cell in each image
Plot_Point_Mapping_tif(lomag_subset, reverse_calc_p_himag_subset, lomag_coord_df, himag_coord_df, lomag_tiles_dir, himag_tiles_dir, _vmax=4000)
#retrieve local coordinates of cells with reads
df_genotype = df_all_reads[(df_all_reads['cell_barcode_count_0'] > 0) & (df_all_reads["well"]==well_ID)].reset_index(drop=True)
lomag_local_coords = np.array(df_genotype[['tile','i_nucleus','j_nucleus']]).astype(int)
#calculate global coordinate of each cell
calculated_P_lomag = Local_to_Global(lomag_local_coords, arranged_grid_lomag, [h_lomag, w_lomag])
#use mapping parameters validated above to calculate global higher mag coordinates of each cell
calculated_P_himag = model_TRS(calculated_P_lomag, mapping_parameters, angle='degree')
#use global higher mag coordinates to translate to local higher mag coordinates, which are associated with phenotype information
reverse_calc_p_himag = Global_to_Local(calculated_P_himag, arranged_grid_himag, [h_himag, w_himag])
#concatenate genotype and phenotype information
df_highmag_coords = pd.DataFrame(data=reverse_calc_p_himag, columns=['tile_himag','i_himag','j_himag'])
df_genotype[['tile_himag','i_himag','j_himag']] = df_highmag_coords[['tile_himag','i_himag','j_himag']]
df_genotype.to_csv("Well"+str(well_ID)+'_reads_himag-mapped.csv')
df_genotype.sample(10)
#the higher mag coordinates can now be used to print out images of cells all unified by some feature described in df_genotype, e.g., sgRNA sequence